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Interferometric gravitational wave detectors are devoted to pick up the effect induced on masses 
by gravitational waves. The variations of the length dividing two mirrors is measured through a 
laser interferometric technique. The Brownian motion of the masses related to the interferometer 
room temperature is a limit to the observation of astrophysical signals. It is referred to as thermal 
noise and it affects the sensitivity of both the projected and the future generation interferometers. 
In this paper we investigate the relevance of small non-linear effects and point out their impact on 
the sensitivity curve of interferometric gravitational wave detectors (e.g. VIRGO, LIGO, GEO. . . ) 
through perturbative methods and numerical simulations. 

We find that in the first order approximation the constants characterizing the power spectrum 
density (PSD) are renormalized but it retains its typical shape. 

This is due to the fact that the involved Feynman diagrams are of tadpole type. 

Higher order approximations are required to give rise to up-conversion effects. 

This result is predicted by the perturbative approach and is in agreement with the numerical re- 
sults obtained by studying the system's non-linear response by numerically simulating its dynamics. 

PACS numbers: 02.50.-r, 05.40.+j 



I. INTRODUCTION 

It is known that every correlation function of a stationary stochastic process may be fully obtained as the functional 
derivative of the associated generator. 

This one has to obey some equations connected to the equations of motion which govern the dynamics of the 
system. If we consider a stochastic system and apply to it all the mathematical techniques and strategies, which 
have been currently set up for studying a particle physics problem, then we get some new ways of approching the 
Fokker-Planck equation connected to a random process, both through perturbative methods and non perturbative 
ones. This requires a procedure of quantization starting from a classical Lagrangian associated to the theory itself 
and a passage to the Euclidean metrics. 

The probability of shifting from one place to another in the case of a Brownian motion is expected to obey an 
evolution equation, like a propagator of Feynman R 

In accordance with that analogy the equation of Fokker-Planck is equivalent to a Schrodinger equation in Euclidean 
time and the proper Hamiltonian operator 
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This means that it is possible to study a stochastic process using the methods of quantum mechanics. Let's introduce 
the a la Schwinger functional 



1 For an exhaustive survey see Jl|-[| and references therein 
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Z(J,K)= J <q,tf\UjK{tf,ti)\q i) ti>dq 
Uj K (t f ,U) = Texp / (J(i)g(t) - iK{t)f>(t))&t 



with q(t) and p(t) taken as operators in the Heisenberg's picture. The equations of motion 

are known as Schwinger-Dyson equations H and rarely they are completely solvable. It has been shown in the 
literature that approximate solutions may be yielded by expanding order by order the above equations if the leading 
term in the drift force is linear. Indeed Z{ J, K) has the form of a path-integral with 

L = E( 4 ^(? - - f>>) + E ^ w 

a, 6 a 

the Lagrangian function in the general case of more than one dimension. What is the physical role of Dl It stands 
for the noise's amplitude and if it tends to zero we recover the classical macroscopic equations without fluctuations. 
It is the same situation as in quantum mechanics where the uncertainty is related to the constant Ti. Exploiting this 
analogy we find an equivalence between 
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dtSJ a (t) jaK SJ{t) 



Z(J, K) = [2D abJ — + K a (t)]Z(J, K) 

W^ hJa(<) + ^ b/b( ^) )]z(J '^ 



and the Heisenberg equations. 

Starting from the lowest order approximation one may solve the Schwinger-Dyson equations above with an iterative 
method. The next step is using the functional 

Z(J,K)=expW(J,K) 

to recover from W(J, K) the connected correlation functions. On the other hand it may be expected that the same 
averages may be obtained following a Langevin approach that is studying a stochastic differential equation whose 
solution is distributed with P(q, t\qo, to) satisfying the Fokker-Planck equation. It may be assumed that after a 
transient period the system becomes stationary and its statistical properties are independent from the initial values 
qo and to- 
ll. MATHEMATICAL FORMALISM 

A physical system is usually represented by many deegres of freedom with their characteristic evolution times. 
Those variables that evolve very rapidly may be integrated over and their contribution is averaged over long time 
periods. What we are interested in is the resulting evolution of the system. The influence of the microscopic deegres 
of freedom on the macroscopic ones may be expressed as a term of interaction |^],[| . 

The characteristic times of such interaction are very short in comparison with the characteristic times of the slowly 
evolving deegres of freedom and the two time scales may be treated separately, at the expense of including a stochastic 
term in the equation of motion written for the macroscopic variables. Such term must be interpreted as follows. At 
any time there are so many interactions that we may only treat them by their statistical properties and consider the 
effect on the macroscopic system as little fluctuations over the deterministic motion. 

Another effect is that phenomenological constants appearing in the equation of motion are related to the interactions 
with the microscopic deegres of freedom. 

So it is only a matter of convenience separating deterministic from stochastic contributions in the equation Jt]Q . 

As an example we treat a classical harmonic oscillator subjected to stochastic forces. It is convenient to use the 
notation of the Langevin equation that is a common equation of motion with a random contribution. The friction 
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constant and the elastic constant are 2ro£wo and mug, m is the mass and D is the diffusion constant whose relation 
with KbT we do not use yet. 
The evolution equation is 

mq(t) + 2m£uj q(t) + mu)%q(t) = F(t) < F(t)F(t') >= 2DS{t - t') . 

The familiar method for solving equations of that type is by means of a vector 

1 \ 



$ = -r$ + F $ = ( 

p 



with the following notation 



p = mq u R = coo^T^e • r = ( 2 £ q ) 



There are two linearly independent eigenvectors characterized by the condition 

p = -muJolS, ± i )q . 

By the similarity transformation 

1 1 -1 ft I jUR\ 

/l_(^ + ^)2 -™ (5 + !^) 1 



the state vector may be expressed in terms of the eigenvectors of T. The eigenvalues we get are 

Formally the solution of the Langevin equation for given initial conditions is 

$'(t) = e- r,t $'(0) + / e r ' (T ^F'(T)dT (2) 
Jo 

where we have introduced 



$' = a/- 1 * = 1 ( q + + ^ p 



mcj (Z + i^)q+p 



wo 



F' = M-\ " ) = 



o , i / -±-(e+ia&)F(t) 



F(t)> ./rr^^ji V n*) 



as state vector and forcing term in the basis of eigenvectors of T. This approach seems rather artificial but it is 
justified by the problem of writing the evolution operator. As a consequence of the transformation we applied above 
the evolution operator becomes e~ r * and its form is now very simple. 

The solutions of the Langevin equation are merely linear superpositions of the eigenvectors of T. The correlation 
functions for the stochastic term are 

<F(*(t)F;(t')>=2D' i:j S(t-t') (3) 

and it should be pointed out that the diffusion constant changes into a tensor |l2^ ]. 
The explicit form is 

, _ Dujq (ml y 2 ;d^(£ - *73fK 

ij 2lu r { -L-Ot + iOB.) 1 ' 

1X mwo vs cjo J 
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and as a consequence 

< i { ) j{ ' >_ 2& R {muJo) 2{ m^oe-^l*- 4 '! (m^o)^-^*-*') ' ' 

Since we are interested in stationary properties t ad t' are very large. If this were not true the initial conditions should 
be taken into account. Even if we preserve the terms containing $'(0) they all vanish for large times. An intuitive 
idea may be drawn by considering the two components in (^|). 

The term in $'(0) doesn't correlate with anything. When we take an average using the properties (Q) the terms in 
$'(0) are simply a factor. 

When times increase they tend faster to zero than possible averaged terms they are multiplied for. The correlation 
functions written above may be turned again in the original basis 

£> e -£w |t-t'| / C^j \ 

< q(t)q{t') > = 3 cosu R {t -t') + — smLO R \t - t'\ 

JJ e -Z^o\t-t'\ 

< q(t)p(t') > = —— ainuj R {t - t') 

_£) e -£w |t-t'| 

< p{t)q{t') > = — — sinw fl (i - 

2m^U)oU) R 

< p(t)p(t') > = JL e -**>|t-*'l (coBU R (t - t') - — sinwfl|i - t'\) . 

Z^LOq \ U) R J 

When we have a connection of several oscillators the procedure stated above is the simplest one because of the 
dimensionality of T. For now it may seem a tortuous mathematical way. 



A. Generating Functionals 

It may be shown that the results of the preceding subsection coincide with the treatment in the introduction. 
This is an example in which an analytic solution is available. If T is a matrix whose eigenvalues are distinct and 
non zero we may use the convention 

a. i 

with the Hamiltonian operator 

// y^D.MM, • ;«!>:. = -1^ . (4) 

i,3 3 

It may be shown that the related Schwinger-Dyson equations we should solve are 

d , SW(J,K) SW(J,K) 

^ + r ^^jwr = Ki(t) + Ai ^w 

d. _ S SW(J,K) 



from which we derive 



i dt 6ij - Tji) -5K^- = - Ji{t) 



d 6 2 W(J,K) 
{ dt Sij+Tij) 8J j (t)5K h (T) =Sih6{t - T) 

d J 2 iy(J,g) 
( d^ " T ^8K 3 {t)5J h {r) = -^ h6{t - T) 



d_ 5 2 W(J,K) _ 8*W(J,K) 

dt ° ij +Lll) SJ 3 (t)SJ h (r) lUii 5K (t)SJ h (r) 
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and the solutions are 

?2 



8Ji(t)8J s (T) ^^ a + Xp 1 j[ 



K(t-to) e -^{r-t )l 



where we have introduced 

D°0 = Q'^Di^f } (5) 
and and f/ are the initial and final times. Using the above results we may write the functional 

tf S 2 W(J,K) f f f f*t S 2 W(J,K) 
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if the reference system is chosen in such a way that &i(t ) — 0. Otherwise it appears a term linear in J that takes 
the initial conditions into account. 

Hence we have constructed a solution for the connected correlations' generator. 



B. Non-linear Corrections 



Now suppose that the response of the oscillator is not perfectly linear. We may write the differential equation 

mq\t) + 2m£u) q(t) + mu)%q(t) + eq 3 (t) = F(t) 

to describe such a situation. The parameter e is assumed so small that the system's behaviour can be described by 
a perturbative approach. On a formal level we are searching for the eigenvectors of the evolution operator that have 
the property of being non-interacting eigenstates of the system. Because of the perturbative term in the equation of 
motion, some new cross-terms appear in the evolution operator whose representative matrix is no more diagonal, but 
for the description of the system the basis of imperturbed eigenvectors is still valid [p||Tc|] . 
The estimated corrections are 



A < $£(t)$i(0 > 

A<$?(t)^ 2 (0 > 
In the physical basis 

A < q(t)q(t') > = 

A < q(t)p{t') > = 
A < p(t)q(t') > = 
A < P (t)p(t') > = 
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are the correlation functions. 

These relations require no knowledge of t and t' but only of t — t'. If t and t' were not very large this property 
should no longer be true. Throughout we use infinitely large values for times, to calculate correlation functions in 
equilibrium conditions. In the stationary limit we may write the power spectrum density (PSD) as follows 

r+oo 

S(f)= e- 2 " /T <q(t)q(t + r)>dT (6) 



since < q(t)q(t + r) > is time-independent and depends on r only. We have 

10, D 2 j ug 

S (f\ = D l _ UC (m U g)'(2m^o)^2^ fro > , ? s 

~ 2(m^ ) 2 (f - )2 + 4^/2 " [(^ _ 2gi)2 + 4 ^2 /2]2 

which can be written 

e m = ^ 1 fsl 

UJ m 2 K 2 (l + _ 47r 2/ 2 ) 2 + (47T/m^ ) 2 1 j 

in the first order approximation. Now we have several ways of interpreting this result. We have searched for sta- 
tionary solutions of the Langevin equation when non-linear terms appear. These terms are related to non-quadratic 
contributions in the Lagrangian function [l4[ . 

This means that in the equilibrium state the proper frequency changes with the local form of the potential function. 
Consequently the parameters governing the responsive behaviour of the system appear modified. As outlined above the 
term in e must contribute very little in the motion equation. Without going into details we say that mui 2 q >> eg 3 . 
Since in equilibrium conditions < q 2 >^ 2£m 2 u 3 ^ e es thriated provision 2 ^mfr^ <<; 1 must be satisfied. This 
parameter is the same appearing in the corrections recovered in this section. They should be wrong if 2 £mfr i « 1 

were not provided. If the system is in its transient state < q 2 > depends on its initial conditions and its value may 

2 

be larger than the constant — -. 

In such a situation we cannot use perturbative methods to study the non-linear component of the force. 



C. Renormalized Coefficients 



Suppose that the oscillator is in contact with a bath and that the probability distribution in the asymptotic limit 
is related to the Boltzmann statistics. We shall show that the dynamic correlation functions we obtained and the 
averages we are getting in the next calculations are consistent. In the stationary limit 

2 J^J^e-^^-Q^^qMqdp KrT 3e K B T 
<q >— 5 3 = 5-tl — 7 -ttt) t — > 00 

J — 00 J- 00 H F 

r r e-T^ { ^ + -^ W+ ^ ) p 2 dqdp 
. 2 ^ J — 00 J— 00 f 1 f 

<p >— 5 j = mKB-L t — > 00 

in the first order approximation. If we call D = 2m^ujoK bT we obtain the same values estimated above for 
< q(t)q(t + t) > and < p(t)p(t + r) > if r = 0. 
We may conclude that 

^ 2^ S- 00 S- o e ~ 1 ^ T ^ d P ^ 2^ JZcf Zo e ~^P 2d Q d P 
<q >= — — — — H <P >= — — — — h 

f°° f°° e K B T dqdp f°° f e K B T dqdp 

are the same with 

rr P* 1 2 2, ? 4 rr P*. 1 2, 1x ,^,2 

= 2m 2™** + £ T ° r H= 2^ + 2 WWo(1 + 3& ^4 )q 
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in the first order approximation. Finally the relation D — Ira^uo^K qT is a thermodynamic property which we expect 
by the fluctuation-dissipation theorem jyj. 

Wherever there is damping there must be fluctuations. They are small because of the factor KbT. Such a 
size is known beforehand because the fluctuations are related to ordinary equilibrium statistical mechanics. From a 
macroscopic point of view the fluctuation-dissipation theorem tells us that the equilibrium is a sort of balance between 
two opposing tendencies (the noise and the damping) and that the fluctuations are determined by the temperature. 
This fact is expressed by the identity D — Im^ujQK bT for the macroscopic parameters. F(t) is treated as an external 
force whose stochastic properties are given by both the constants 2m(uo and T and it acts regardless of q. Even 
if F(t) is irregular and unpredictable we may fully express its properties in terms of the postulated autocorrelation 
function. Then the differential equation 

mq(t) + 2m£_u} q{t) + mv%q(t) = F(t) < F(t)F(t') >= 2D6(t - t') 

is well-defined with a well-defined solution. Its characteristic Green function is obtained from inserting 



F(t) = S(t) and has the form Q(t) = M e -C"o*^«* 



(9) 



where we intended u>r — luq y/l — £ 2 as before. Let us introduce 

mq{t) + 2m&oq(t) + muilq(t) + eq 3 (t) = F(t) 

This Langevin equation defines 



/-hoc r+oc 
g(t - t)F{t)At -e Q{t- r)g 3 (r)dr 
-oo J — oo 

/+oo /"+oo r p-\-oo />+oo 

Q(t-T)F(T)dT-e 9(t-r) / G(t - r')F(r')dr' - e / Q{r - rV(r')dr' 
-oo J — oo LJ — oo J — oo 
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/-foo r-\-oo r r-\-oo 

g{t-T)F{ T )dr-e G(t — t) / g{T-T l )F{T l )dT l 

-oo J — oo Vy — oo 

/+oo /"+00 /"+oo 

Q{t-t')[ g(r' - r")F{r")dT" - e / S(r' - r") 9 3 (r")dr"] 3 d7 
-oo •/ — oo «/ — oo 



At 



At 



which is clearly a series in e as it should. For fixed small values of e we may consider the first terms only. The 
initial point is q(t) if e = 0. It has an autocorrelation function < q(t)q(t') >= R(\t — t'\). It is permissible to use the 
expansion in e in order to get the corrections 

/+oo r+oo 
G(t - T)R(\t' - T \)R(0)dT - 3e / Q(j> - r)R(\t - r|)i?(0)dr 
-OO J — OO 

/+oo r+oo 
g(t - r)R(\t' - t\) / g(r - t')R(\t - t'\)R(\t - T '\)R(0)dTdT' 
-oo J —oo 
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/-roo r-roo 
g(t' - T)R(\t - r|) / g{r - t')R(\t - t'\)R(\t - r'|) J R(0)drdr' 
-oo J —oo 

/+oo /*+oc 
g(t- T )R(0) / g(r-r')R{0)R(\t' -r'\)dr'dT 
-OO J —OO 

/ + 00 /•+00 
g(t'-r)R(0) / £(r-r')i?(0)i?(|i-r'|)dr'dr 
- oo ./ — oo 

/+oo /-+oo 
£(t-r) / £(r-r').R(|r-r'|) 2 .R(|i'-r'|)dr'dr 
-oo J —oo 

/+oo r+oo 
S(t' - r) / 0(t - t')R(\t - T'\) 2 R(\t - r'Ddr'dr 
-oo J —oo 

/+oo /"+00 
0(t-r) / g(t' -t')R(0) 2 R(\t- r'Ddr'dr 
-oo — oo 

/+oo /-+oo 
S(t-r) / £(£'-r').R(|r-r'|) 3 dr'dr 
-oo J —oo 
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III. NON-LINEAR EFFECTS 



After the preliminaries of the previous section we transfer our attention to a single pendulum and moreover we use 
the second order approximation to estimate the corrections related to a non-linear term in the motion equation. Thus 
if we use same notations and techniques as before, we find the power spectrum density 

4£tJg K B T 

+ — ~ 9a > - U \ + (A^OWj [i- ~ —(j) I - 2(1-^) Ip^p^+t^ + (^-9^) 2 + (6?a^) 2 J 

with uj = 2tt f and a = J Kb T\1 . From the Boltzmann distribution in the non- linear case we easily recover the result 

< n >= J-ocJ-oc = („_!)!!. AflJ , ,« 



J- co J- ooe -l ? iT(^ + ^9 2 +^) d(zd p K) n/2 2 2(1 +2H)* J 

This relation is a sort of scaling-law in the second order approximation. As an example we may write 



< 9 > — o < g >< 3 >= — 3a(— — jt) + 14a ( -^— 2 j 



'(tow 2 ) 2tow 
and compare it with 

< q\t) > conn = -3a(^) 2 + 14a 2 (^) 2 . 
This is the asymtpotic expansion of < q A (t) > CO nn for i — > oo. In the same notation 

r.+00 

e v 2 

— OO 



,K B T j^e (V + t?' 4 )^ 2 ^' KbT 3a 



<q >= n-) ^ = o- 1 - + 6a ) 



It is also possible to obtain a numerical evaluation of the involved integrals and make a comparison of the result 
with the expansion in a. One has for q 2 two series of points as in Fig.[| for many choices of a. 

We put q' = y j^j^ g in the integrals. It's necessary for the expansion to be valid that the value of a be small. For 
a great it's logical that the perturbative approximation is not longer justified. 
We also know that 



< ? 2 (*) >= S(f)df = -^+6a' 



from the corresponding dynamic correlation function recovered by (|^) . 

We are less interested in static properties than in dynamic properties, so turn our attention to the power spectrum 
density. We may apply the previous analytic estimations to a simple pendulum having length I ~ lm and mass 
to ~ 500kg. 

Indeed even in a multi stage pendulum holding the mirrors of an interferometer, the thermal noise depends on 
the last pendulum and on the test mass itself, and above the pendulum proper frequency the dynamics mimics the 
response of a free mass. In order to have an idea of the noise if some non-linearity arises we use a value for to staying 
between a single test mass and all the tilting filters in the VIRGO superattenuator Jl3| . 

Accordingly to the antennas' pendula requirements that are being built, we take I as the coarse order of magnitude 
for the length of the system's suspension. 

This argument inspired the choice of the mechanical quality factor's value too. 

In Figjl] the PSD is shown. It gives the analysed system's behaviour. 

Because of the non-linear contributions, there are two corrections particularly peaked in correspondence of 2itf ~ f 
and27r/-3f. 

Another consequence of non-linearity is related to the bins correlation and we develop it in this section. Owing to 
white noise if the system were linear bins would be completely uncorrelated. On the contrary from the non-linear 
components of the response it follows that 

< i<?(/)i w)i 2 > conn = -^df - 27r2(/ e iy ,2) )|5(/)|2|5(/,)|2 ^ • 
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Here 2r\ is the width of each bin. Thanks to the smallness of KbT it can be found that the bins correlation is 
dramatically small. We report the related plot in Fig.|| with the usual coefficients as before. If both / and /' 
approach the resonance frequency the bins correlation shows its top value. 



IV. DIAGRAMS REPRESENTATION 



As a companion to a perturbative calculation of the solution of the stochastic differential equation, it is sometimes 
advantageous to consider the alternative method of stochastic diagrams. 

The solution of the linear differential equation is associated with a line and a cross denoting the Green's function 
and the stochastic force respectively. If the non-linear contribution is taken into account we may iteratively solve the 
stochastic differential equation. The non-linear term contributes a vertex that represents the convolution product of 
the unperturbed Green's function and the solutions of the equation. 

This procedure is better illustrated by the first diagrams in Fig.^. 

The notation is meant to indicate the following. From the differential equation we may write the solution as an 
iterative expansion. Accordingly q(t) depends on various products of noise terms which we average over to obtain 
correlation functions. When there are several equivalent possible ways of combining crosses a combinatorical factor 
arises. 

The first order correction to the two-point correlation function is the result of two tadpole-type diagrams. Their 
contribution is a correction to ujq. 

In the second order approximation there are tadpole-type contributions as well. 

Yet there are two- loops diagrams that introduce new terms with different poles than the previous corrections. 

Now observe the diagrams in Fig.|. The first three ones contribute a second order correction to loq. The other two 
entail new terms which are peaked near the resonance frequency and three times the resonance frequency. If £ is large 
the former peak may not be distinguished from the principal one. 

The latter is evident. 

The elastic constant becomes 

3 

mt^Q — > towq(1 + — ol — 9a 2 ) 

and the damping constant 

27a 2 

2m£w -> 2to£w (1 - y^j) 

In order to understand those modifications we may note that there are stronger forces pulling the system back to the 
equilibrium configuration and the result is that the measured recall constant towq has a new greater value. 

For linear systems in contact with a heat-bath it seems that the damping force have been expressly designed to 
constrain every motion below the thermal noise. 

This will be true even in the non-linear case but due to the system's dynamics and the appearance of up-converted 
proper frequencies damping times are longer. 

This effect participates to the attenuation of the system's response below the fundamental frequency as it is shown 
in Fig.Q. The coupling constant assumes non-perturbative values and the numerical results confirm the tendency of 
the perturbative case. 



V. NUMERICAL SIMULATIONS 



In this section we are planning a numerical simulation to analyse the equation 



q + 2& q + Lo 2 q = 2J KbT(u} ° F(t) (10) 



with D = 2<m& K B T and < F(t)F(t') >= 6(t - t') as previously stated. 
For convenience we write 



p = -2£iv p — moj^q + V2DF \ p 
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= + < EifiiEWj >= 6(t - t 1 ) ( 





4m£GJ K B T 



where /i(<&) is not derived from a potential. This may be immediately seen from the cross derivative difj ^ djfc. The 
stochastic term S is completely characterized by its zero mean and correlation matrix. This latter one is singular. 

If we look at the original equation of motion we see that its solution depends both on the initial conditions and the 
distribution of the noise source F(t). 

As a consequence the statistical properties of the system are characterized by 

P(q",p",t"\q',p',t') =< q",p"\U(t",t')\ q ',p' > 
with \q' ,p' > simultaneous eigenvector of q and p. The evolution operator is such that 

U(t", t') = e -^(*"-*') H = DU 2 p - ifl p (muj 2 q + 2& p) + ^fl q p 

and the following rules 

[ii p ,p] = -i [ii q , q] = -i [n p , q ] = [ti q ,p\ = [n p , ry = [p, q ] = o 

are to be satisfied. 

In other words q and p are regarded to as independent variables. At any t their value is distributed according to 
the probability that a certain initial configuration has evolved into another one because of a stochastic force whose 
statistical properties are known. 

We may also recover the distribution of every physical quantities depending on q and p. At every time the 
corresponding operator may be expressed by means of the Heisenberg picture 

Cl(q(t")p(t")) = U- 1 {t\t')h{q{t'),p{t'))U{t" 1 t') ±{l( q (t)p(t)) = [H,n(q(t),P(t))} 

if the observable quantity fi does not explicitly depend on t. Making a comparison with the Poisson brackets of 
classical mechanics we note a reversed sign. Utilizing the Baker-Campbell-Hausdorff formula in U(t",t') one obtains 

P(q",p",t"\q',p',t') = J dgi f dq 2 J dpi j dp 2 < \q 2 ,p 2 > x 

< q2,P2\e-( t "- t ' ) " 2 \qi, P2 >< qi^e'^^W ,p' > +0 ((*" - t'f) 
with H = H i + H 2 . If it happens that raw' 2 = we have 

1 (v" - v ' e -2S"o(t"-t'))2 

P(p",t"\p / ,t / )= - =exp [-^ p e L- ] (11) 

because of the identity 

cxp[-(t" - t')(DUl - i2& flpp)} = exp[--^- (1 - cxp[-£w (*" - t')}) ft?] x exp[*2£wo(t" - t')tl p p] 

Next for the purpose of simplifying the algebra we write t" — t' = At. Now we examine the situation specified by 
mio 2 7^ 0. Finally we derive a three-steps process to O(At) 2 order. Successive higher orders can also be obtained. 
We first demonstrate that 

5(q"-q') ( (p" -P' + mLU 2 q'At) 



<q", P "\e-^ Hi \q',p' > = ^= ;; cx P i- 



V2DirAt V 4£,A * 
Hi = DTI 2 — imoj 2 Tl pq 

< q",p"\e- At ^\ q ',p> > = 5(q"-q' - 5-^— [1 - exp(-2^ Ai)]) 5{p" - p' exp(-2^ At)) 

H 2 = -i2£ujoTlpP H Tl q p 

m 



starting from showing that 



10 



<q",p"\e-^ Al \q',p' > = <</'.;/ , 
< q",p"\e~^\q\p' > = < q",p"\q'(At),p'(At) > 



^\q',p' >=<q",p"\e 



f 


#2,9 


= -P 

rn ^ 




^2,P 


= -2& p 



<?(t) = ±p'(t) q'(0) = q' 

p'(t) = -2& p'(t) p'(0)=p' 



- DTT 2 i / / 2 



where it may be muj^q — > miJ^q + eg 3 . We yield the following algorithm 



<7i = <?' 



p 2 = {l-2£,uj a + \{2^ Mf)p 1 

92 = 9i + 2^(2£a;oAi - ±(2^ Ai) 2 ) 



P"=P2 + #/(<? 2 ) + V^Al6 
q" = Q2 



where /(g) = muj^q + eg 3 but any /(g) could be substituted. 

With this method an algorithm of any precision may be planned. The fluctuations of the physical quantities are 
very small due to D = 2m£uioK bT '. Because of this we shall scale both the physical variables p — > y/KsTmp and 
1 ~~* \/ ^rmJ 2 ~ Q ^° S e t two variables with mean square value one. 



VI. RESULTS AND CONCLUSIONS 



In the following we shall use the same £ and luq previously adopted. 

These are the only two important physical parameters because of the scaling we introduced at the end of the 
previous section. The form we recovered for S(f) is to be multiplied by the scaling factor ■ It is convenient to 
recall that we wrote it as the spectrum at the zero order with a term in the denominator standing for the "self-energy" . 
It gives a correction to the physical parameters. 

Furthcmore there is a brand new resonance appearing as an up-conversion of the original single resonant frequency. 

In FigJ^ the power spectrum density for a = 0.1 is shown. There is a little peak at three times the resonance 
frequency as previously expected. 

The results obtained from the numerical simulation are compared with the curve referred to a = 0.0. It may be 
noted that S(f) is reduced for low values of / and that the fundamental peak moves a little forward. Indeed we 
expected this from the perturbative variations of muj^ and 2m^ujo. 

The up-converted resonance is set in a frequency region corresponding to small values of S(f). This is due to the 
response of the system that behaves like a filter attenuating those components whose frequency is greater than ojq. 

For smaller values of a the up-conversion effect is not appreciable. Let us discuss if a — ^^Jy, << 1- Actually 

the quotient n^jrp is generally tremendously smaller than the parameter (2KbT)^ 1 . The conclusion is that the 
coupling constant between non-linearities in the suspension system's behaviour and thermal noise is so small that 
an enhancement of up-conversion effects is not realistic. This is due to the smallness of thermal fluctuations and 
to the fact that there are not contributions proportional to the coupling constant, because only tad-poles diagrams 
contribute at the first order approximation. These diagrams contribute a parameters' correction but do not modify 
the correlations at different times. 

We expect that the resonant frequency should be greater if a goes up. It happens if the system explores the 
non-linear zone of the phase space. When T rises the greater amplitude of fluctuations just produce such an effect. If 
the recall constant increases it is obvious that the opposite tendency is obtained. 

Finally we note that formally there are two more peaks than in the linear case, but one of these can not be 
distinguished from the fundamental one whose width is proportional to £. This is related to nothing but the PSD's 
proper shape and its typical parameters. 

APPENDIX A: EVOLUTION OPERATOR DECOMPOSITION 

In this section we want to give a brief proof of the exactness of the identity 

exp[-(t" - t')(DYl 2 p - i2^ ft P p)] = exp[-— — (1 - exp[-fwo(t" - 01) x exp[*2£wo(t" - t')%p] 
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First of all we multiply the evolution operator by the identity such as to get 

exp[-(i" - t')(Dtil - i2& %p)] = exp[-(i" - t')(DIL 2 p - i2& fl p p)} exp[^— n^] exp[-^— 
and than we use 



D -21 r D ~ 2 . 



(t"-t')Hfr2(t"-t')H _ fr2 



to write the form 

e -(t"-t')H e ^n 2 Pe {t''-t')H e -(t''-t')H e - 1 ^ 5 ti 



n^exp(-^ (r-t')) 



s_fi 2 — 



with H = DTI 2 — i2^cJoii p p. The commutation relation [7J,n p ] is sufficient to have the equation (Al). 
It simply describes the temporal evolution in the Heisenberg picture. The final identity is yielded using 

DU 2 . r DU 2 . 



inside the (A2). Substituting the expression for H we arrive at the following simple identity 



(Al) 



(A2) 



that is the outcome. 

There are many other ways of demonstrating the above identity. Here we cite the iterative way of developping the 
Baker-Campbell-Hausdorff formula in E^j. 
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FIG. 1. Power density spectrum for a pendulum with a quality factor ~ 100 




FIG. 2. Bins correlation < \x(f) 2 \x(f') 2 > fo „. The width of each bin is 277 = O.lmHz and the adimensional coefficient is 
a = 0.1. Such coefficient represents the level of non-linearity of the concerned system 
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FIG. 3. The expansion in a leads to good approximations of < x 2 >. The approximated expression is not longer valid when 
a — > 1. There is an asymptotic limit < x 2 >~ ol^iKbT^lo 2 ,) -1 if a grows 
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FIG. 4. Graphical representation of x(t) and < x(t)x(t') >. When the average over the random force is taken, all crosses 
are joined together in all possible ways. For example one can get < x(t)x(t') > by combining two crosses at a point in the 
corresponding product 
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FIG. 5. Second order corrections for the two-point correlation function. There are contributions including tadpoles. A 
tadpole is nothing but a contraction of x with itself at the same time. Since < x(t)x(t) > is a constant tadpoles are insertions 
of constant terms 
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FIG. 6. In this figure the imperturbed theoretical curve related to the linear case and the power spectrum density corre- 
sponding to a numerical simulation of the non-linear one are both shown. Data are considered over a period T chosen as 
~ 1200^-. In the frequency domain the variations due to non-linearities are concentrated in a band that is as narrow as Q is 
high. It seems that for Q — > oo even small non-linearities may be revealed 
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FIG. 7. In this figure a sketch of how the sensitivity curve for a pendulum is modified by a non-linear term is shown. As an 
example an effective reduction is achieved at low frequencies. The non-linear force causes the pendulum to oscillate with two 
important modes. In order to keep small the amplitude of the up-converted fundamental frequency a (which is proportonal to 
T) is expected to be less than 0.1. These results are in accordance with the estimation of our perturbative calculations 
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